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ABSTRACT 

Aims. We explore the lensing properties of asymmetric matter density distributions in Bekenstein's Tensor- Vector-Scalar theory 
(TeVeS). 

Methods. Using an iterative Fourier-based solver for the resulting non-linear scalar field equation, we numerically calculate the total 
gravitational potential and derive the corresponding TeVeS lensing maps. 

Results. Considering variations on rather small scales, we show that the lensing properties significantly depend on the lens's extent 
along the line of sight. Furthermore, all simulated TeVeS convergence maps strongly track the dominant baryonic components, non- 
linear effects, being capable of counteracting this trend, turn out to be very small. Setting up a toy model for the cluster merger 
1E0657 - 558, we infer that TeVeS cannot explain observations without assuming an additional dark mass component in both cluster 
centers, which is in accordance with previous work. 
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1. Introduction 

As is known, General Relativity (GR) cannot explain the dy- 
namics of our universe on large physical scales as the amount of 
visible mass clearly lies below what would be expected from 
the observed gravitational effects. Commonly, this is denoted 
as the Missing Mass Problem. The usual remedy is to invoke 
a form of matter which does not couple to light, therefore be- 
ing referred to as Dark Matter (DM). Over the last years, this 
paradigm has been remarkably successful in forming a consis- 
tent cosmological picture since it provides suitable explanations 
for the observations on supernovae la ([Astier et al. 2006 1, large- 
scale structure (Tegmark et al.|[2004{|Cole et al. |2005] >rand the 
CMB HSpergel et al.||2007) . 

However, one can also take a different point of view 
and modify the law of gravity itself. In the past, there have 
been several suggestions for such modifications: f(R) gravity 
( |Starobinskij| 198Q]>, c onformal Weyl gravity and Aether-type 



theories (IZlosnik et al. 



alization of the latter, t 



2007 1, to name just a few. A special re- 



ity (TeVeS) (Bekenstein 2004 Zlosnik et al. 



te so-called Tensor- Vector-Scalar grav- 
" 2006l l has re- 



cently gained interest as it provides a fully relativistic frame- 
work for the Modified Newtonian Dynamics (MOND) paradigm 
(Milg rom| |1983| |Bekenstein & M ilgrom 1984). Compared to 
other modifications, MONDian dynamics is characterized by an 
acceleration scale oq, and its departure from classical Newtonian 
predictions depends on acceleration: 



fi \ — la = -V<b N + S. 



(1) 



Here, <S>n denotes the common Newtonian potential of the visible 
(baryonic) matter and S is a solenoidal vector field determined 
by the condition that a can be expressed as the gradient of a 
scalar potential. The function fi, controlling the modification of 



Newton's law, has the following asymptotic behavior: 



fi(x) 
fi(x) 



X X <K 1, 
1 X » 1. 



(2) 



Eq. ([TJ has been constructed to agree with the fact that the ro- 
tation curves of spiral galaxies become flat outside their cen- 
tral parts (Sanders & McGaugh, 2002). Analyzing observational 
data, Milgrom estimated ao ~ 1 x 10~ 10 m/s 2 . Within this work, 
we shall study TeVeS and its built-in MONDian dynamics in the 
context of gravitational lensing, focusing on non-spherical den- 
sity distributions. 

The paper is structured as follows: Starting with a brief in- 
troduction to TeVeS and the formalism of gravitational lensing, 
we shall have a look at analytic lens models in spherical sym- 
metry. In TeVeS, the effect of gravity on matter is controlled by 
a function which is only constrained by its non-relativistic lim- 
its, i.e. the Newtonian and the MONDian limit, with the func- 
tion's actual form in the intermediate part providing an addi- 
tional degree of freedom to the theory. Therefore, we shall pro- 
ceed with a general analysis of this so-called free function's in- 
fluence on the TeVeS deflection angle. Choosing a specific form 
of the free function, we will present a numerical tool that al- 
lows the treatment of non-spherical lenses in TeVeS. In contrast 



to already existent MOND solvers (Ciotti et aL| 2006 Brada & 



Milgrom 1995] |199"9| our method is based on fast Fourier tech 
niques, achieving high-resolution solutions for the TeVeS scalar 
potential on time scales up to a few hours on standard PCs. 
Highlighting basic properties and limitations of this method, we 
shall apply it to a set of different matter density distributions in- 
cluding discussions on the failure of the thin lens approximation 
and general properties of the TeVeS lensing maps. Finally, we 
will create a toy model for the cluster merger 1E0657 - 558, 
with the resulting lensing maps confirming and extending the 
basic findings of Angus et al. (2007]!. 
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2. Fundamentals of TeVeS 



2.2. The Free Function 



In the following, we will give a brief review on TeVeS and the 
approximations used for quasi-static systems like galaxies (or 
galaxy clusters) and cosmology. If not specified in any other 
way, the "square" of a vector denotes the square of its Euclidean 
norm, i.e. (A) 2 = (||A|| 2 ) 2 . Henceforth, we shall use units with 
c= 1. 



2. 1 . Fields and Actions 

TeVeS gravity is based on three dynamical fields: an Einstein 
metric g^ v , a vector field such that 



g» v U,U v = -1 



(3) 



and a scalar field <p. An essential feature of TeVeS is the intro- 
duction of a physical frame described by the metric which is 
needed for gravity-matter coupling only and obtained from the 
non-conformal relation 

lw = V - 2U M U y sinh(2<£). (4) 

The geometrical part of the action is exactly the same as in GR: 

1 



16ttG 



(5) 



where is the Ricci tensor of g^ and g the determinant of g MV . 
Note that the TeVeS constant G must not be mistaken for the 
Newtonian gravitational constant Gn (cf. Sec. 2.3 1. The vector 
field's action S v reads as follows: 



K 



32ttG 



F^F, 



A(g" v U^U v 



(6) 



with F^y = U^ v — U YJ1 . Here the constant K describes the vector's 
coupling to gravity and A is a Lagrangian multiplier enforcing 
the normalization given by Eq. ([3]). Eq. (|6]l corresponds to the 
classical Maxwell action, the field now having an effective 
mass. The action S s of the scalar field (f> involves an additional 
non-dynamical scalar field cr, and takes the form 



c 4 

cr 2 F», v + -^FikGcr 2 ) 



(7) 



where W v 



U^W and F is a dimensionless free func- 



tion. As the field cr is related to the invariant ht" <pji(p, v , however, 
it could in principle be eliminated from the action. While k is 
the coupling constant of <p to gravity, the constant I is related to 
Milgrom's ao and has the dimension of a length (see Sec. 2.2 1. 
Finally, according to the equivalence principle, the matter action 
is given by 



get X. 



(8) 



Matter fields are coupled to gravity by the physical metric g^ v , 
i.e. world lines are geodesies of the metric g^ rather than g MV . 
As usual, the corresponding equations of motion can be derived 
by varying the total action S = S g + S v +S s +S m w.r.t. the basic 
fields. 

In order to obtain Newton's law in the non-relativistic high 
acceleration regime (a » qq), the coupling constants k and K 
have to be small, i.e. 



(9) 



Therefore, TeVeS is kept close to GR in a sense that it will re- 
cover well-known features of GR, albeit modified by the other 
fields. 



In TeVeS, the transition from Newtonian dynamics to MOND is 
controlled by the free function F. Following Bekenstein (2004), 
the "equation of motion" for the non-dynamical field cr suggests 
introducing a new function p.(y) which is implicitly given by 

1 



• (*F(M) ~ ~H F (ji) = y, 



(10) 



with 

kGa 2 =^kl z hr<l> 4l <P tY )=n(y). (11) 

For further analysis, we shall assume the function fi(y) to behave 
well in a physical sense, i.e to be smooth and monotonic in both 
cosmological (y < 0) and quasi-static situations (y > 0). In or- 
der to reproduce both a MONDian and a Newtonian limit, the 
quasi-static branch of the inverse function y{fj) has to satisfy the 
following conditions: 



y(jf) -> oo 
yQi) ~ bfi 2 



H -> 1, 
fi «c 1, 



(12) 



where b is a positive real constant. If this is the case, the constant 
I can be related to Milgrom's oo by 

<h~k 

°° = AnZl * 4^' (B) 
where S = 1 - K/2 - 2<f> c and <p c is the cosmological value of 
the scalar field which is assumed to be small (<p c <k 1). In Sec. 
4. 1 we shall return to the free function and its properties in the 
context of gravitational lensing, concentrating on the branch rel- 
evant for quasi-static systems. 

2.3. Quasi-static Systems 

According to |Bekenstein| ( |2004 1, the physical metric field near a 
quasi-static galaxy (cluster) is identical to the metric obtained in 
GR if the non-relativistic gravitational potential is replaced by 



(14) 



where <& N is the Newtonian potential generated by the baryonic 
matter density p. In this approximation, it is consistent to take 
JJ^ = (U' , 0, 0, 0) which can be shown from the corresponding 
field equations. Then we have 

kfhT^y -> kl 2 (V4>) 2 (15) 

and the equation of the scalar field reduces to 

V [/j (kl 2 (Vcf>) 2 ) V0] = kGp. (16) 

Eq. ( [T6] > corresponds to the non-linear elliptic boundary value 
problem and can be treated numerically. In Sec. 5.2 we shall 



give a detailed description of the method we use to determine 
the solution for the scalar field <p, including a discussion on its 
problems and limitations. 

Since we have K, <p c <k 1, the quantity S has a value close to 
unity, and the total potential <5 can essentially be written as the 
sum of the common Newtonian potential and the additional 
scalar field, i.e. Eq. ( fT4| ) may further be reduced to 

<D = <Div + 0. (17) 

As the constant G is related to the Newtonian gravitational con- 
stant Gn by (Bekenstein] [2004} 



- + 4^) G ' 



(18) 



we will additionally assume G ~ Gn throughout this work. 
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2.4. Cosmology 

Similar to the case of GR, it is possible to derive a cosmological 
model in TeVeS. Assuming the basic fields to partake of the sym- 
metries of the Friedmann-Robertson- Walker (FRW) spacetime, 
the analog of Friedmann's equation reads 



— ipe *+ p j- — + - 
3 a 1 3 



where is the energy density of the scalar field given by 



^ ^ V by •> 



-2/iyQi)+/i 2 FQi) 
4k 2 PG 



(19) 



(20) 



Since we are interested in the physical metric, we have to make 
use of transformation (EJ and finally obtain 



1 da j, (a ■ \ 
a at \a I 



(21) 



with 



dt = e^dt, 



= e^a. (22) 

In order to simplify matters, however, we shall introduce the 
"minimal" cosmological model proposed by Zhao et al. (2006): 



According to |BekensteTn| ( |2004| l, it is consistent to assume 
that the cosmological scalar field evolves slowly in time through- 
out cosmological history. Thus, its contribution to the Hubble 
expansion is negligibly small, with a ratio 0{k) compared to the 
matter contribution. Setting = and recalling that <p <K 1, the 
physical Hubble parameter can be expressed as 



H 2 *H 2 * H 2 (o m (l + z) 1 + Q A + Q K (l + z) 2 ) , 



(23) 



where D.k ~ 1 - £2a - £l m . Since there is no DM in TeVeS, 
we have to consider a minimal-matter cosmology that should be 
consistent with observational data in order to obtain a reason- 
able cosmological model. |Zhao et al.] actually find a good fit of 
the high-z SNe distance moduli data set by choosing an open 
cosmology with Q A ~ 0.46, Q. m ~ 0.04 and H ~ 70km/s/Mpc. 
However, they also point out that when moving to very high red- 
shifts, this open cosmology has problems, i.e. it underestimates 
the last scattering sound horizon, which actually seems to be an 
artifact of the crude approximation as recent work has shown 
(Zhao 2006 1. Nevertheless, in the context of gravitational lens- 
ing, this simple model is sufficient for assigning the distances of 
lenses and sources up to a redshift of z ~ 3. 



3. Gravitational Lensing in TeVeS 

3. 1 . Light Bending in Slightly Curved Spacetime 

In general, light rays move along the null geodesies of the under- 
lying metric field, i.e. the null geodesies of the physical metric 
g^y considering the framework of TeVeS gravity. For an on av- 
erage homogeneous and isotropic universe with local perturba- 
tions, however, light mostly travels through unperturbed space- 
time and is only deflected close to inhomogeneities which act 
as gravitational lenses. If the non-relativistic potential O and the 
peculiar velocity v of the lens are small (<L>, v <K 1), we can pre- 
sume a locally flat spacetime being disturbed by the potential 
these conditions are well satisfied for galaxies and galaxy clus- 
ters. 

It is well-known ( Schnei der et al.| |1992| l that the deflection 
angle of a light ray under these assumptions can be expressed as 



- 2 J" 



V ± <M/ = a CR + 2 V^dl, 



f 

•J — cx 



(24) 



where <t> is given by Eq. (fl7j, V ± denotes the two-dimensional 
gradient operator perpendicular to light propagation and integra- 
tion is performed along the unperturbed light path (Born's ap- 
proximation). In addition to the deflection angle caused by the 
Newtonian potential Ojy, there is a contribution arising from the 
scalar field <j>. Because <p is connected to the matter density in a 
highly non-linear way, it is not possible to relate the projected 
matter density to a two-dimensional scalar deflection potential 
just like in GR (cf. Sec. 3.2 1. Therefore, we are obliged to solve 
Eq. ( fT6l > for calculating the TeVeS deflection angle, which is a 
very delicate issue (cf. Sec. 5.2 1. Compared to the distances be- 
tween lens and source and observer and source, however, we 
may still assume that most of the bending occurs within a small 
range around the lens. This enables us to fully adopt the GR lens- 
ing formalism which will be introduced in the next section. For 
further discussion, we shall choose coordinates such that unper- 
turbed light rays propagate parallel to the z-axis. 



3.2. Lensing Formalism and Critical Curves 

In gravitational lensing, it is convenient to introduce the deflec- 
tion potential *P(0) ( |Schneider et al.||1992) : 



V(0) 



D s D d J 



®(D d 0,z)dz, 



(25) 



where we have used 6 — ^/D d . Here £ is the 2-dimensional 
position vector in the lens plane, and D s , D d , and D ds are the 
(angular diameter) distances between source and observer, lens 
and observer, and lens and source, respectively. If a source is 
much smaller than the angular scale on which the lens proper- 
ties change, the lens mapping can locally be linearized. Thus, 
the distortion of an image can be described by the Jacobian ma- 
trix 

BP _ (l - K - J\ -J2 

-y 2 1 



80 



k + 7i 



(26) 



where fi=i]/D s and 17 denotes the 2-dimensional position of the 
source. The convergence k is directly related to the deflection 
potential *P through 

k = i VP (27) 



and the shear components j\ , 72 are given by 



1 

71 = 2 



(d Zx ¥ 

ae 2 



d 2 ^ 

do 2 , 



72 = 



d 2x V 



2 + r 2 2 - 



Points in the lens plane where 
det^l = (1 -K- 



y)(l - K + y) = 0, 



(28) 



(29) 



form closed curves, the critical curves. Their image curves lo- 
cated in the source plane are called caustics. Since images 
near critical curves can significantly be magnified and distorted, 
which, for instance, is indicated by the giant luminous arcs 
formed from source galaxies near caustics, a thorough analy- 
sis of the behavior of such curves in a TeVeS universe will be 
profitable. 

3.3. Analytic Solutions 



Following Zha o et al7| ( |2006| l, we switch to a notation which turns 
out to be more suitable for analytic studies. Instead of the func- 
tion ix, we shall consider a new function ft which is defined by 



4n 



I-P 



(30) 
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Fig. 1. The TeVeS Hernquist lens (dashed) with r# = 2kpc 
compared to its corresponding point lens (solid) w.r.t. the de- 
flection angle a for M = 10 n M o , a = 1 x 10~ 10 m/s 2 and 
D = D ds D d /D s = 850Mpc. For f » 1, both angles approach 
the constant = 2tt s/GMao m 0.58 . The transition to the 
MONDian regime can be characterized by the critical radius 
r = y/GM/ao ~ lOkpc. 



where k, K are the coupling constants of the scalar field <f> and the 
vector field U^, respectively. Similarly, we can relate the func- 
tion y to another function 6$ in the following way: 



4n 



(-!) 



(V0) 2 



(31) 



where b is the real-valued parameter of the function y(p) in ( p~2| >. 
Choosing the free function such that 



In the limit r# — > 0, the Hernquist lens coincides with a point 
lens. In this case, we find that the deflection angle is given by 



4GM , 

&(£) = — ^— + In VGMa . 



(36) 



Obviously, the scalar part of Te VeS gravity of a point mass seems 
to mimic the presence of a dark isothermal sphere. Therefore, 
both GR including DM and TeVeS will essentially make the 
same lensing predictions for £ being much larger than the exten- 
sion of the lens, but the highly non-linear coupling of the scalar 
field strongly suggests that there may be significant differences 
when moving to the strong acceleration regime near the center. 
Note that, although Eqs. p4| and ( [36]l do not explicitly depend 
on k and K, a is still given by Eq. ( |13[ ). 

Fig. [T] shows the lensing properties of both the TeVeS 
Hernquist lens (r# = 2kpc) and its corresponding point lens 
where we have set M = 10 n M o , a = yjbk/4nl = 1 x l(T 10 m/s 2 
and D = Dd s D c i/D s = 850Mpc. Since its deflection angle can be 
expressed analytically, the Hernquist lens is a perfectly suitable 
candidate for testing an algorithm for non-spherical problems. 

4. Influence of the Free Function 

Considering a spherically symmetric situation and applying 
Gauss's theorem for a spherical surface of arbitrary radius, Eq. 
(fToTl can be transformed into 



= - — VOa,. 
4-n/i 



(37) 



Assuming we already know cDai, for example by solving 
Poisson's equation, the relation above can directly be used to 
calculate V</> for any given function fi(y) (Remember that y = 
kl 2 (\7<p) 2 for quasi-static systems, with I given by Eq. (\3\). If ft 
or (f> N cannot be obtained analytically, treatment with numerical 
methods, which can easily be applied in the spherically symmet- 
ric case, becomes necessary. Because of their simplicity, spheri- 
cally symmetric systems are particularly suitable for investigat- 
ing the effects of the free function y(fj) on the deflection angle. 



d-/i) 2 



(32) 



it is possible to obtain an analytic expression for the deflection 
angle of a Hernquist lens, i.e. a lens whose matter distribution 
follows a Hernquist profile (Hernquist 1990| l given by 



Mr, 



2nr(r + r H y 



(33) 



the Hernquist radius r# being the core scale length and M the 
total mass. Eq. ( [33] > is a spherical profile which closely approxi- 
mates the de Vaucouleurs 7?^ 4 -law for elliptical galaxies. Using 
elementary calculus, we eventually end up with 



where 



r H A(£) 



4£ VGMao ■ 



4GMg 



4GM£ 



(34) 



arsinh ^|l - (r H /^f\ % < r H 
arcsin yi - (r H /£) 2 £ > r H 



(35) 



4.1. Parameterization of the Free Function 



Having set the cosmological background in Sec. 2.4 we shall 
focus on the free function's quasi-static branch (y > 0): If y(fi) 
can be analytically continued into the ring domain R = {z e C | 
< \z — 1| < 1}, it can be expanded into a Laurent series. Thus, 
y(jl) takes the following form for < fi < 1 : 



n=l V w n=0 



(38) 



with coefficients a n , b„ € R. Expanding the above expression for 
fi <K 1 to second order, we must have the following relations for 
the coefficients a„, b„ to keep the second condition in ( fT2| ): 



b 2 + 



b + ^ a„ = 0, 

oo 

b\ + ^ a„n = 0, 

n=l 

n(n + 1) 



(39) 



±0. 
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Fig. 2. Absolute (left) and relative (right) diffe rence between the TeVeS deflection angles of the analytic and the numerical Hernquist 
lens choosing the parameters from Sec. 



3.3 



Clearly, the deviations are small, Aa <; (10 4 ) and A&/a ana i v < (10 4 ). For the 
numerical calculation, we have assumed k = 0.01 and y(fi) = n 2 1(1 - jj). 



As a simple example, we take the function 

,- 2 



and find that the non-zero coefficients are given by 
a\ = 1, bo — —I, b\ — — 1. 



(40) 



(41) 



Setting the coefficients a„ and b„, we are able to directly control 
the specific transition behavior from MONDian to Newtonian 
dynamics. Using the expansion ( |3~8j ), we shall study the effects 
of varying the free function y(ji) on the deflection angle within 
numerical analysis. 

4.2. Comparison to the Analytic Model 



Taking the simple choice ( |40) , we compare the numerical result 
of the TeVeS Hernquist lens to the analytic solution derived in 
Sec. 3.3 Fig. 2] shows the absolute and relative difference be- 
tween the deflection angles of the analytic and the numerical 
lens models using the same parameters as in Sec. |3.3| where we 
have set k = 0.01 for the numerical calculation. Although we 
have considered two different free functions y(fi), the actual dif- 
ferences are fairly small, Aa < (10" 4 )" and A&/a ana i y < (10~ 4 ), 
and thus negligible with respect to today's observational resolu- 
tion limit. 

Considering the above, it appears that we could be able to 
determine classes of free functions y(p) that nearly produce the 
same deflection angle. Although we do expect an infinite num- 
ber of such classes, a closer look will be of advantage (cf. Sec. 
|5.1) . For a systematic approach, we shall make use of the pa- 
rameterization introduced in Sec. 4. 1 to analyze the effects of a 
varying free function. 

4.3. Varying Parameters 

Still considering deflection by the Hernquist lens, we now focus 
on a varying form of the free function y(jx) and the value of the 
coupling constant k. For our investigation, all remaining param- 
eters are chosen as in Sec.|3.3|unless specified in any other way. 



Starting from an arbitrary set (a n ,b„), e.g. the set pT) , we 
begin with a variation of the constants a„ where we have to adjust 
bo, b\ and bj according to ([39]), setting k to a fixed value, e.g. k = 
0.01. Interestingly, numerical analysis has shown that there are 
no significant changes to the deflection angle for a wide range of 
parameters, i.e. the relative changes are comparable to those we 
found in the last section, A& re i < 10~ 5 - 10~ 4 . For instance, if we 
choose 

aig = 1, ai9 = -2, a20=l, (42) 



which is the expansion of 



(1-/0 



20' 



(43) 



the relative deviation is of (9(10~ 5 ). As we have seen, the sin- 
gularity at fi = 1 recovers the Newtonian limit in quasi-static 
situations, but remarkably, the transition behavior seems almost 
insensitive to the particular realization of this singularity. 

Similarly, we have also examined the effect of a chang- 
ing coupling constant k taking the coefficients a„, b„ as con- 
stants. Again, the relative differences turned out to be very small, 
Aa re i < 10~ 5 - 10~ 4 , varying k within the range of 10~ 4 - 10~ 2 
for different sets (a„,b n ). Obviously, as long as it is small, i.e. 
k < 0.01, the calculation of the deflection angle does not really 
depend on the exact value of k. 

As for the coefficients b„ with n > 2, however, there is a 
strong influence on the deflection angle, basically allowing to 
create arbitrary transitions from MOND to Newtonian dynamics. 
In accordance with the above analysis, it seems that the b n alone 
can be used to characterize the free function. In general, the exact 
form of y(f£) has to by constrained by observational data being 
independent of the particular law of gravity, which is subject to 
other work, e.g. Zhao & Famaey (2 006| l. 



5. Non-spherical Lens Models 

Within this section, we will examine the properties of more gen- 
eral lens systems using numerical methods. Introducing our al- 
gorithm for the treatment of non-spherical lenses in TeVeS, we 
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will investigate a set of different matter distributions including a 
toy model of the cluster merger 1E0657 - 558. 

5. 1 . Choice of the Free Function 

Setting k = 0.01, we shall restrict all further analysis to the fol- 
lowing form of y(ju): 



d-A') 2 ' 



(44) 



We will make use of this specific y(p) for two reasons: First of 
all, the choice (j44|» is easily inverted, i.e. 



(45) 



and therefore it is possible to express the derivative with respect 
to y analytically: 



1 



dp 

dy = 2^(1 + y/y) 



2 ' 



As will become clear in Sec. |5.2| both p(y) and dp/dy are part of 
Eq. ( fTo} . Since dp/dy — > oo fory — > 0, a possible solver of (jT6]» 
might be extremely sensitive to the corresponding run of dp/dy 
in that regime. By choosing the analytic expressions ( |4"5] l and 
{46) , respectively, we are able to avoid numerical inversion and 
differentiation of the free function, which may prevent a desta- 
bilizing influence on the algorithm. 

Secondly, our choice allows us to use the analytic Hernquist 
lens for comparison in order to test the accuracy of a non- 
spherical solver for this specific density profile. According to 
Sec. [4] Eq. ( [44] > is close to the choice ([32} and produces nearly 
the same deflection angle, thus justifying such a comparison. 

5.2. Calculating the Scalar Potential 

Since p = p(y) and y = kl 2 \V<p\ 2 for quasi-static systems, an 
expansion of the l.h.s. of Eq. ( fT6) i yields 

2^kl 2 ((dMdjMdidrfj) + = kGp. (47) 



Defining an effective matter density p such that 



where 



kG kl 2 dp i \ 
P=—p-2—£ ((drfXdjMdidjcf))) , 



p dy 



(48) 



(49) 



we may choose an appropriate first guess of (/> and calculate an 
initial density p (0) by using Eq. ( |49) . Solving Poisson's equation 
by means of Fourier methods, i.e. Eq. ( |4"8"j ) with the r.h.s. being 
fixed (p = p (0) ), we find a new field which can be used to 
obtain p (1) and so forth. 

Without any further modification, this approach fails to con- 
verge in most cases, with the tjf n) oscillating rapidly. Including 
a relaxation into the iteration, however, it is possible to enforce 
convergence for a variety of problems, and thus our final iterative 
scheme reads as (p (0) is calculated from an initial guess <po) 



A^ (n) = p (n \ 
t (n+r> = tj4> (n) + (l -u)cf) in \ 



(50) 



where we have introduced the relaxation parameter co e 
additional iteration field and 



p-2 



(«) ; „;2 



(»)_ kG „ J d n kl 



dy) p ( 

(«) 



an 



(51) 



^=n(y {n \ (|) =|(y (B) ). y (n) = kiW\ 

For suitable values of co, our method turns out to work very well 
for a wide range of density profiles (cf. Sec. 5.2.3 1. However, 
our investigation has shown that the relaxation's success is very 
sensitive to the particular choice of w, i.e. u> has to be chosen 
from a very narrow range, co = 0.75 + 0.1. Although conver- 
gence is achieved within a wider range of co, its behavior quickly 
deteriorates. Fortunately, this value seems to be almost indepen- 
dent of the particular density profile, and therefore it will not be 
necessary to adjust to once it has been determined for a certain 
density. 



(46) 5.2.1. Point Lens Approximation 



As the scalar field's gradient decreases much more slowly com- 
pared to the Newtonian one far away from the lens, one would 
actually be obliged to move to very large volumes in order to 
neglect contributions from outside the box and obtain correct re- 
sults for the deflection angle. Thus, assuming a fixed grid size, 
this would excessively degrade the resolution of the correspond- 
ing two-dimensional lensing maps. In the following, we shall 
discuss an approximation allowing us to avoid this problem: 

Considering a finite grid with N + I points per dimension 
(jV is chosen as an even number), we may rewrite the scalar part 
of the deflection angle as the sum of contributions coming from 
both inside and outside the grid's volume: 

XfA.v ^oo 
V ± <f> m dz + 4 V x <t> (out) dz, (52) 
fA.v JfA* 

with the quantity Ax denoting the distance between neighboring 
grid points. Assuming that the scalar field at the boundaries is 
approximately given by that of a point lens, i.e. 



dout) 



VGMao log(r), 



(53) 



we obtain the following expression (M denotes the total mass 
inside the volume): 



f fA, 

a s = 2 V ± <; 



(in) 



dz + 4A, 



where 



VGM 



(NAx 

arc tan 

I 2q 



(54) 



(55) 



and q 2 = x 2 + y 2 . If applicable, we need to perform the integra- 
tion only over our finite grid since all contributions from outside 
the box can be expressed analytically. Concerning our iterative 
solver, we may additionally assume the boundary conditions of 
the fields to be of spherical symmetry, and it turns out to be 
sufficient to use ( |53| l as an initial guess for (p. The highest reso- 
lution compatible with the limitations of our computer hardware 
is set by N - 384 which is used for all numerical calculations. 
Before turning to non-spherical lens systems, however, we shall 
examine our method's accuracy. 
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Fig. 3. Accuracy test of our numerical method: Assuming a Hernquist lens with the parameters from Sec. 
numerical result for the deflection angle to the analytic one. The left panel shows the relative difference Aa re \ = (&analy - 
&num)/&anaiy Please note that the quantity Aa re \ is limited by -0.02 < Aa re i < 0.01 for reasons of presentation, values outside 
this range are truncated. The right panel illustrates the relative deviation for the central part where Aa re i reaches a maximum of 
approximately 20%. 



5.2.2. Accuracy 

Comparing the numerically obtained deflection angle of a 
Hernquist lens to the analytic result ( |34"| ), we will determine the 
accuracy of our tool assuming the parameters from Sec. |3.3| 
which corresponds to a galaxy-sized mass distribution. As pre- 
viously mentioned, such a comparison is justified according to 
our analysis in Sec. 4.1 Concerning the numerical setup, we 
choose a grid volume of V — (50kpc) 3 (the lens is placed in 
the grid's center), and in order to obtain a sufficiently large value 
of D = Dd s Dd/D s , we set the redshifts of source and lens to 
Zsource = 3 and ziens = 0.63, respectively. If not explicitly noted, 
we shall keep this choice of redshifts throughout the following 
sections. At £ « 25kpc, the relative deviation between the ana- 
lytic Hernquist and its corresponding point lens is approximately 
6%. Although this difference is quite large, we assume the va- 
lidity of Eq. ( |53) ) in order to balance accuracy and resolution 
(N = 384). Thus, we are able to resolve structures up to a mini- 
mum extent of Ax « \3Qpc. 

The left panel of Fig.[3]shows the relative difference Aa re \ = 
(flanah - &num) I &anah between the numerical and the analytic de- 
flection angle of the Hernquist lens. For reasons of presentation, 
we have limited the range of Aa re i to -0.02 < Aa re i < 0.01. 
Ignoring the very center of the map (right panel), we find the rel- 
ative deviations in the interior are of (9(10~ 3 ). Moving outwards, 
i.e. to larger £ = |£|, these deviations increase and reach values 
up to 5 - 6% at the grid's boundaries (£ > 25kpc). However, as 
long as £ < 15kpc, we still have Aa re i < 1%, again neglecting 
the central part. The large differences close to the boundaries are 
likely to be a mixture of artifacts caused by the Fourier transform 
of actual non-periodic fields and contributions due to ( |55] l which 
become more significant with increasing ^. Having a look at the 
right panel of Fig. [3] we see that Aa re i strongly increases in the 
central region reaching a maximum value of roughly 20%. The 
reasons for these large deviations are probably related to both 
the limited resolution of our grid and the small values of a in the 
center (see also Sec. 5.2.3 I. According to Sec. 3.3 the TeVeS 
deflection angle of the analytic Hernquist lens decreases to zero 
for f — > 0. Since this transition happens on a rather small scale, 



our numerical model cannot fully recover the deflection angle in 
the central region. Moreover, we have to consider that the matter 
density ( |33] l becomes infinite at £ = 0, which, of course, cannot 
be accomplished in a numerical calculation. Due to the grid, this 
singularity is smoothed out, causing an effective loss of mass in 
our numerical model. This loss has an overall influence on the 
deflection angle and may significantly contribute to the errors we 
have discussed above. Investigating non-spherical systems, how- 
ever, we shall only consider lenses which follow smooth density 
distributions. 

For numerical simulations that similarly allow using the 
point lens approximation, we may assume an accuracy equal to 
that of the Hernquist lens. As we are mainly interested in the 
strong lensing regime, we shall restrict ourselves to the grid's 
interior where the relative deviations are of<9(10~ 3 ). Due to fi- 
nite resolution, however, we expect the accuracy to degrade to 
some extent in regions where the deflection angle a approaches 
values close to zero when moving to more generic lens systems. 
Although smooth density profiles will probably not produce de- 
viations as large as we have found around ^ = for the Hernquist 
lens, we cannot make any specific statements on the quality of 
our simulations in such areas. Still, this should affect but a frac- 
tion of the overall result, thus being acceptable for the following 
analysis. 



5.2.3. Problems 

As for the solver of the scalar field, we have encountered the 
following problems: Considering more complicated density dis- 
tributions, we have found the relaxed iteration to be less effi- 
cient, i.e. the iteration generally takes more time to converge. 
As it turns out, this cannot be compensated by changing the re- 
laxation parameter at, which would actually lead to even worse 
convergence properties or a complete failure of the method. Still, 
the additional amount of time that has to be employed is accept- 
able in most cases and corresponds to a factor of 2 - 3. 

Furthermore, independently of the particular value of co, we 
encounter the relaxation to generally fail for certain choices of p. 
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Fig. 4. Left: Numerically calculated Te VeS convergence map for the King-like profile ( |56] ) expressed in terms of the GR convergence 
K gr assuming to - 50kpc. Right: Since ( |56] l corresponds to an axisymmetric configuration, the effective TeVeS convergence k and 
the ratio K/K gr depend on the radial coordinate £ only. As the calculated convergence maps relatively deviate from circular symmetry 
by <9(1(T 3 ), which is due to our Fourier method, the presented results are averaged over all directions. Choosing zo = 50kpc (solid) 
and zo = 400kpc (dashed), we see that k is significantly amplified in the central region when moving to higher values of zo- 



Since it can mostly be resolved by slightly modifying the orig- 
inal density profile, this second problem is probably of purely 
numerical origin. However, we point out that it might also hint 
on an exceptional behavior of the scalar field <p that is not ac- 
cessible to our solver. Considering the Hernquist profile, for 
instance, the difficulties found in the central part may reflect its 
intrinsic instability w.r.t. TeVeS/MOND, rather than a negative 
feature of our code ( jCiotti & Binney|[2004| . 

5.3. Thin Lens Approximation 

As our first task, we want to investigate the validity of the thin 
lens approximation in TeVeS. According to former work consid- 
ering lensing in classical MOND (Mort lock & Turner] |2001| l, we 
expect a break-down of the approximation due to the non-linear 
coupling of the scalar field to the three-dimensional matter den- 
sity. In the following, however, we are rather interested in quan- 
tifying this break-down by exploring the lensing properties of a 
mass distribution being contracted or stretched along the line of 
sight, i.e. the z-direction, making use of our new numerical tool. 
For this purpose, let us consider a three-dimensional density dis- 
tribution p following a King profile ( King 1972 1 which is given 
by 



p(f) = Po 



1 + 



2\ 



(56) 



where r c is the core radius and po me matter density at r — 0. Eq. 
|56) is an empirical law that fairly describes the distribution of 
both galaxies and gas inside a galaxy cluster. However, in order 
to analyze TeVeS effects which are only due to the lens's extent 
along the line of sight, we have to parameterize its "thickness" 
and additionally ensure a constant projected mass density. Thus, 
we introduce a slightly modified profile: 

p(q,Z)=p Q(q)Z(z), (57) 

where 



Q(q) = 



1 + 1-1 

<7o 



(58) 
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Fig. 5. Radii of the inner (dashed) and outer (solid) critical curve 
for different choices of zo: Increasing the lens's extent along the 
line of sight, the radius of the inner (outer) critical curve is driven 
outwards showing relative changes of up to roughly 6% (16%). 
The critical lines are calculated by interpolation between the grid 
points. 



and 



Z(z) 



1 + 



2\ 



(59) 



with qo, zo > being the corresponding core lengths and 
q 2 = x 2 + y 2 . Since the expressions ( |58] > and |59) are obtained 
by integrating ( |56]l o ver one and two dimensions, respectively, 
our new choice ([57]) is actually kept close to the original King 
profile. Varying the parameter zo, we are now able to directly 
control the lens's extent in the z-direction. 



M. Feix et al.: Asymmetric Gravitational Lenses in TeVeS and Application to the Bullet Cluster 



9 




\x [kpc] \ x [kpc] 



Fig. 6. Numerically calculated TeVeS convergence k (left) and the corresponding ratio K/K gr (right) for an elliptical profile of type 
( |60| l with Xo = 350kpc and yo = Zo = 200kpc: In the central part, the TeVeS convergence k is unevenly amplified, thus breaking the 
elliptical symmetry. Compared to the dominant GR contributions, however, this effect is almost negligible. 



Concerning the numerical setup, we take the point of origin, 
i.e. x — y — z = 0, to coincide with the grid's center and as- 
sume the density profile ( |57| > to be smoothly cut at a radius of 
R - 1.5Mpc, which is necessary as, otherwise, our King-like 
profile would contain an infinite amount of mass. In accordance 
with the point lens approximation, we set the grid's volume to 
V = (5Mpc) 3 , thus having a spatial resolution of Ax m 13kpc 
(N = 384). Furthermore, we take g = 200kpc and choose po 
such that the lens's total mass is given by M = 10 15 M Q , which 
now corresponds to a cluster-sized mass distribution. Regarding 
the remaining parameters, we keep the settings introduced in 
Sec. - 



3.3 



the redshifts of source and lens being fixed to z SOU rce - 3 
and ziens = 0.63, respectively. Henceforth, if not explicitly noted, 
all presented numerical results are based upon these settings, in- 
cluding the spherical cut-off of the particular density profile at 
R = 1.5Mpc. As we have to meet condition ( f53j ) in order to 
apply the point lens approximation, we are obliged to satisfy 
Zo ^ 600kpc in our numerical simulations. 

The right panel of Fig. |4] shows the effective TeVeS conver- 
gence of our King-like profile expressed in terms of the corre- 
sponding GR convergence K gr for zo = 50kpc and zo = 400kpc, 
respectively. Note that the GR maps are independent of the par- 
ticular choice of zo- Increasing the value of zo, we observe a 
significant amplification of the TeVeS convergence around the 
center while there is basically no change in the outer region. As 
expected, the variation of zo has no effect on the symmetry prop- 
erties of the convergence map. Concerning the TeVeS shear map, 
we find a similar behavior: While there is a strong increase of y 
in the very center, we find only small changes in the outer parts. 
Interestingly, the TeVeS shear is not exactly circularly symmet- 
ric in that region any longer, with the actual form depending on 
the particular extent of the lens. Rather than being intrinsic to 
TeVeS, however, this is probably due to Fourier artifacts caused 
by the scalar field solver or the point lens approximation, an in- 
fluence of the latter being actually expected as the choice of zo 
has an impact on Eq. ( |53"j ). 

Let us continue our analysis considering the effects on the 
critical lines due to the changes of k and y: Since |57) is axially 
symmetric, the corresponding lines turn into circles. In Fig. [3] 
the radii of both the inner and outer critical curve are presented 
for different values of the parameter zo- Obviously, these radii 



are increased when stretching the lens along the line of sight, 
showing relative deviations of up to roughly 6% and 16% for 
the radial and tangential critical radius, respectively. Note that 
the critical lines are calculated by interpolation between the grid 
points, thus allowing to determine positions which are below the 
grid's resolution. Although our investigation is limited to a small 
range of zo, we find appreciable differences between the lensing 
maps which are assumed to considerably grow when stretching 
the lens further. 

Summarizing the above, we may conclude that the lens's ex- 
tent along the line of sight significantly affects the strong lensing 
properties. Therefore, the mass distribution along the z-axis can 
be regarded as an additional degree of freedom in TeVeS. 

5.4. Elliptical Lenses 

In this section, we shall consider lens systems whose projected 
mass density follows an elliptic profile. Therefore, introducing 
the scale lengths xo,yo,Zo > 0, let us consider a matter density 
distribution of the form 



P(r) = Po 



1 + 



xo 



2\ 



(60) 



Keeping yo and z fixed, yo - Zo - 200kpc, we investigate the 
lensing properties for different choices of the parameter xo, again 
setting the total mass to M — 1O 15 M . Fig. [6] illustrates both 
the TeVeS convergence k and the corresponding ratio k/k„, with 
Xo set to a value of 350kpc. Although the symmetry properties 
of the GR convergence map are virtually sustained in TeVeS, 
we can observe an interesting feature located in the central part: 
Compared to its neighborhood, there is a slightly increased am- 
plification close to the semi-major axis, breaking the elliptical 
symmetry. If this effect was larger, it could actually account 
for loosing track of the baryonic matter distribution, thus yield- 
ing a qualitatively different looking TeVeS convergence k. For 
xq = 350kpc, we additionally present a simulation where the 
density profile |60} has been rotated around the z-, y- and x- 
axis by 10°, 20° and 30°, respectively. Clearly, the ratio K/K gr 
illustrated in Fig. [7] (left panel) shows essentially the same in- 
ner structure as in Fig. [6] Therefore, it seems unlikely that the 
observed effect is a numerical artifact caused by our method. 
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Fig. 7. Left panel: Numerically calculated TeVeS convergence ratio K/K gr for the rotated profile ( |o*0| ) with xo = 350kpc: In the central 
part, the TeVeS convergence k is unevenly amplified, similar to the unrotated case shown in the bottom panel of Fig. [6] Thus, it 
seems unlikely that the observed effect is a numerical artifact caused by our method. Right panel: Critical curves for both TeVeS 
(solid) and GR (dashed) assuming an elliptical profile of type d60b with xq = 150kpc and yo = zo = 200kpc. 



In the right panel of Fig. [7] we compare the critical lines in 
TeVeS to those obtained in GR assuming xq = 150. The found 
symmetry-breaking effect does not appear to have any signifi- 
cant influence on the critical curves which therefore do not show 
any unfamiliar shapes compared to elliptical GR lenses. As the 
TeVeS convergence k is calculated by a weighted amplification 
of K gr , however, the critical curves appear at a larger distance 
from the origin and their forms are varied to some extent com- 
pared to GR. Varying the value of xq from lOOkpc to 400kpc, we 
substantially obtain the same findings. 

5.5. Lenses with Multiple Components 

Next, we want to explore gravitational lensing by multiple ob- 
jects. For this purpose, let us consider a rather simple case and 
start with two density distributions, p\ and p%, following the 
King profile (|56]>. Choosing r c = 200kpc and Mi + M 2 = M = 
10 15 M o (Mj denotes the total mass of the object located at r, in- 
side our volume), we shall place our densities at the following 
positions inside the grid volume (r = |r| = corresponds to the 
grid's origin): 



(61) 



Thus, varying the parameters x 2 and Z2, we are able to control 
the relative alignment of our objects along the line of sight, i.e. 
the z-direction, as well as perpendicular to it. 



5.5.1. Equal Masses 

As a first approach, we shall assume the total mass M to be 
evenly distributed on our two bullet-like objects, i.e. M\ = Mi- 
Varying %2 from lOOkpc to 300kpc, we calculate the convergence 
maps and critical lines for different alignments along the line of 
sight, with the results for Z2 = 400kpc (and Z2 — for the critical 
curves) presented in Fig. [8] (top and middle panel). Again, we 
notice that k is amplified such that the symmetry properties of 





(x 2 ) 
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r\ = 





, r 2 = - 
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the surface density K gr , are virtually conserved, similar to the re- 
sult found in the last section. Having a look at the central region, 
we additionally observe that k is increased between the object's 
positions, which is actually expected since the Newtonian gradi- 
ent VO/v becomes small there. Altogether, as the TeVeS conver- 
gence map closely tracks the baryonic matter distribution, and 
we do not encounter any new surprising TeVeS effects consider- 
ing our two-bullet system. 

Increasing the quantity Z2, we discover a significant growth 
of k around the central part, which is in accordance with our pre- 
vious result from Sec. |5.3| Consequently, the corresponding crit- 
ical lines, shown on the r.h.s of Fig. 8](top and middle panel), are 
spatially driven outwards. Please also note that, due to the non- 
spherical symmetry of our problem, the shape of those curves is 
slightly changed when varying Z2- 



5.5.2. Different Masses 
In analogy to Sec 



5.5.1 



we can perform a similar simulation 
choosing 3M\ = M2. Assuming X2 = lOOkpc, both the calcu- 
lated convergence map and the critical lines are presented in the 
bottom panel of Fig.[8]for Z2 = 400kpc (and Z2 = for the criti- 
cal curves). As can be seen from the ratio K/K gr , the convergence 
is more strongly amplified in the £ x > regime, i.e. the region of 
lower mass density. Accordingly, the corresponding critical lines 
are drawn further outwards in that region. As the MONDian in- 
fluence increases for smaller values of the Newtonian gradient's 
modulus I VOjv|, however, this is exactly what one would expect. 
Choosing other bullet alignments or mass weightings, we basi- 
cally obtain the same results. 

5.6. Modeling the Bullet Cluster 

Only recently, the cluster merger 1E0657 - 558, has been an- 
nounced as a direct empirical proof of the existence of DM 



(Clowe et al. 2006, Bradac et al. 2006 1 as the weak lensing 



reconstruction of k shows peaks that are clearly detached from 
the dominant baryonic components, i.e. the plasma clouds. 
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Fig. 8. Lensing properties of our two-bullet system assuming X2 = lOOkpc and Mi = M2 (top panel), xi = 300kpc and M\ = M2 
(middle panel), and xi = lOOkpc and 3M\ = M2 (bottom panel), respectively: We present the numerical results for both the TeVeS 
convergence k (left) and the corresponding ratio K/K gr (middle) setting Z2 = 400kpc. On the right, the TeVeS critical lines are plotted 
for Z2 - (dashed) and 12 - 400kpc (solid). Note that the radial critical curve for xq = 300kpc and Z2 - does not appear due to 
the grid's finite resolution. 



Using an analytic model, Angus et al. ( 20Q7\ have fit this 
map and derived the corresponding baryonic matter density in 
MOND-like gravity, concluding that it is not possible to model 
the merger without assuming an additional invisible mass com- 
ponent located in the central parts of the two clusters. As they 
have used the weak lensing reconstruction of Clowe et al. (2006), 
however, their convergence map does not account for the ob- 
served strong lensing features within this system. 

In addition, the same convergence map appears to be incom- 
patible with the values of plasma masses independently esti- 
mated from Chandra observations, which, as remarked by the 
authors, may be due to the smoothing scale of the weak lensing 
reconstruction. To avoid such discrepancies, we shall create a 
simple toy model of the bullet cluster's baryonic matter density 
that allows to be treated with our numerical method. 

The bullet cluster basically consists of four objects: The main 
cluster, a slightly smaller subcluster and two plasma clouds ap- 



pearing in between. For the three-dimensional matter density, 
we shall model these components using the analytic King pro- 
file ( |56] >, choosing core radii r c of 200kpc and 150kpc for the 
clusters and the plasma clouds, respectively. Concerning the par- 
ticular masses, we use the values of |Clowe et al.| ( |20 06) which 
are derived independently of gravitational lensing. According to 
the authors, the plasma mass is reconstructed from a multicom- 
ponent three-dimensional cluster model fit to the Chandra X-ray 
image while, assuming a mass-to-light ratio of M/L t = 2, the 
stellar mass is calculated from the I-band luminosity of all galax- 
ies equal in brightness or fainter than the component's brightest 
cluster galaxies (BCG). Together with the approximate positions 
of the components, the corresponding mass values are presented 
in Tab.[T](|r| = denotes the grid's center). Please note that all 
masses are averaged within an aperture of lOOkpc radius around 
the given position. For the calculation of the lensing maps, we 
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Fig. 9. TeVeS convergence maps for our toy model of the bullet cluster: Assuming the framework of TeVeS, we present the numeri- 
cally obtained convergence k (top panel) and the ratio K/K gr (bottom panel) for the parameter sets A (left) and C (right). 
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Fig. 10. TeVeS maps of the shear components 71,72 f° r our toy model of the bullet cluster: Assuming the framework of TeVeS, we 
present the numerically obtained shear components j\ (top panel) and 72 (bottom panel) for the parameter sets A (left) and C (right). 
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Component 


Position (x,y,z) [kpc]° 


Plasma mass M x (10 12 M o ) 


Stellar mass M, (10 12 M S ) 


Mum (10 12 M o )* 


Main cluster 


(-350,-50,zi) 


5.5 


0.5 


6.0 


Main cluster plasma 


(-140,50,z 2 ) 


6.6 


0.2 


6.8 


Subcluster 


(350, -50,z 3 ) 


2.7 


0.6 


3.3 


Subcluster plasma 


(200,-10,z 4 ) 


5.8 


0.1 


5.9 



" For each component, the position perpendicular to the line of sight is approximately determined from the corresponding Magellan and Chandra 
images. 

Concerning the masses of our toy model components, we use those derived by Clowe et al. 1 2006 1. Note that all masses are averaged within 
an aperture of lOOkpc radius around the given position. 



Table 2. Parameter sets used within the toy model of the cluster 
merger 1E0657 - 558 



Parameter set* 


Zl 


Z2 


Z3 


Z4 


A 














B 





300 





-300 


C 





500 





-500 


D 


300 


100 


-300 


-100 



In our simulations, the above sets are used to specify the compo- 
nent's alignment along the line of sight, i.e. the z-direction. 



set the redshift of the bullet cluster, i.e. the lens, to zum - 0.296 
( |Clowe et al.| |2006 ) and assume a source redshift of z S ource - E 
Since the position of the particular constituents can only be 
constrained perpendicular to the line of sight, we actually have 
substantial freedom in selecting their alignment along the z- 
direction. For our analysis, we choose four different sets of z,- 
which are shown in Tab. [2] Let us briefly discuss the meaning 
of these choices: Clearly, the parameter set A implies that all 
components are located in the same plane perpendicular to the 
line of sight. The choices B and C account for the plasma clouds 
to be displaced in opposite directions along the z-axis, which is 
a reasonable assumption considering today's view of the bullet 
cluster to be a post-merger snapshot. Finally, the parameter set 
D fairly describes the situation of the axis connecting the cluster 
centers being inclined w.r.t. the z-direction. 

Fig.|9]illustrates both the resulting TeVeS convergence map k 
and the corresponding ratio K/K gr for the parameter sets A and C 
listed in Tab. [2] Again, we find that there are TeVeS effects caus- 
ing additional structure within the central part, and increasing the 
constituent's relative displacement along the z-axis, we observe 
these structures growing stronger to some extent. Similar to our 
previously considered lens models, however, the TeVeS effects 
are not large enough to account for displacements from the dom- 
inant baryonic components. Additionally, we present maps of the 
TeVeS shear components j\ and j2 for the sets A and C which 
are shown in Fig. 
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Clearly, our result confirms the findings of | Angus et al.| 
(2007). To provide an acceptable explanation of the observa- 
tions, TeVeS needs an additional mass component centered at the 
cluster positions. As has been suggested by others, e.g. |Sanders| 
(2007), primordial neutrinos with mass on the order of 2eV 
might be able to resolve the problem. Checking the v-values, 
i.e. the arguments of the free function (cf. Sec. |2.2[), near the 



cluster centers in our simulation, we find that the non-linearity 
of Eq. ( fT6] > still has a relevant impact on the resulting scalar field. 
Therefore, it is not possible to isolate neutrino effects as a pure 
additive contribution to the overall convergence map and to give 
constraints on the amount and distribution of such neutrinos for a 
given mass. Current work is trying to find an approximative way 
of dealing with this issue using our previous numerical results. 

6. Conclusions 

In this work, we have analyzed the effects of gravitational lens- 
ing within the framework of TeVeS, focusing on asymmetric sys- 
tems. 

Considering spherically symmetric lenses, we introduced a 
parameterization of the free function y(ji) showing that the par- 
ticular realization of the singularity at = 1 has a trifling influ- 
ence on the deflection angle. Furthermore, we concluded that 
variations of the coupling constant k lead to negligibly small 
effects as long as k < 10~ 2 . Thus, we were able to determine 
classes of y(ji) that nearly produce the same deflection angle. 

Choosing a single form of the free function, we succeeded in 
building a fast Fourier-based solver for scalar potential <p which 
could be applied to a set of different non-spherical lens types. 
Concerning variations on rather small scales, we noticed a strong 
dependence of the lensing properties on the lens's extent along 
the line of sight, with a significant impact on the critical curves. 
Additionally, every simulated TeVeS convergence map showed 
a strong resemblance with the dominant baryonic mass compo- 
nents, other effects, being capable of counteracting this trend, 
turned out to be very small. To study a more complex lens sys- 
tem, we finally created a toy model of the bullet cluster's bary- 
onic matter density. The outcome of our simulation clearly con- 
firms the results of An gus et al.| ( |2007[ ) as it is not possible to 
explain the observed weak lensing map without assuming an ad- 
ditional dark mass component in both cluster centers. 

Future work will address even more complex lens systems 
like, for example, the galaxy cluster Abell 2390 with its straight 
arc. Increasing the lens's level of substructure, effects that ac- 
count for loosing track of the projected matter density could be- 
come more important, thus influencing the ability of TeVeS to 
model such a lens. In addition, one can use these models to check 
if TeVeS is consistent with the assumption of massive neutrinos 
in galaxy clusters. 
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